From nonlocal gap solitary waves to bound states in periodic media 
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Abstract 

o 

Solitary waves in one-dimensional periodic media are discussed employing the nonlinear 
O |. Schrodinger equation with a spatially periodic potential as a model. This equation admits 

two families of gap solitons that bifurcate from the edges of Bloch bands in the linear wave 
spectrum. These fundamental solitons may be positioned only at specific locations relative 
fC) ' to the potential; otherwise, they become nonlocal owing to the presence of growing tails of 

exponentially-small amplitude with respect to the wave peak amplitude. Here, by matching the 
tails of such nonlocal solitary waves, higher-order locally confined gap solitons, or bound states, 
are constructed. Details are worked out for bound states comprising two nonlocal solitary waves 
in the presence of a sinusoidal potential. A countable set of bound-state families, characterized 
by the separation distance of the two solitary waves, is found, and each family features three 
distinct solution branches that bifurcate near Bloch-band edges at small, but finite, amplitude. 
Power curves associated with these solution branches are computed asymptotically for large 
solitary-wave separation, and the theoretical predictions are consistent with numerical results. 

> 

^ ■ 1 Introduction 

VD ' 

Nonlinear wave phenomena in periodic media are currently attracting a great deal of research inter- 
Q\ ■ est in nonlinear optics, Bose-Einstein condensation and applied mathematics (see Christodoulides 

et al. 2003; Kivshar & Agrawal 2003; Morsch & Oberthaler 2006; Skorobogatiy & Yang 2009; Yang 
2010 for reviews). Apart from scientific curiosity, this research activity is also driven by various 
potential applications, ranging from light routing in lattice networks (Christodoulides & Eugenieva 
2001) to image transmission through nonlinear media (Yang et al. 2011). 

A characteristic feature of periodic media is the existence of bands in the linear spectrum 
where linear disturbances, the so-called Bloch modes, may propagate. Between these Bloch bands 
are bandgaps in which linear disturbances are evanescent but nonlinear localized modes, commonly 
known as gap solitons, are possible. Since the first theoretical prediction (Christodoulides & Joseph 
1988) and experimental observation (Eisenberg et al. 1998) of fundamental solitons in the semi- 
infinite bandgap of one-dimensional (ID) periodic waveguides, various types of gap solitons in 
one- and multi-dimensions have been reported theoretically and experimentally. Examples include 
2D fundamental gap solitons, vortex solitons, dipole solitons, reduced-symmetry solitons, vortex- 
array solitons, truncated-Bloch-wave solitons and arbitrary-shape gap solitons (see Yang 2010 for a 
review). Multi-soliton bound states in periodic media have also been constructed in the framework 
of a discrete nonlinear Schrodinger (NLS) model (Kevrekidis et al. 2001). In addition, the stability 
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of some of these gap solitons has been examined (Pelinovsky et al. 2004, 2005; Shi et al. 2008; 
Yang 2010; Hwang et al. 2011). 

The plethora of gap solitons in periodic media calls for systematic identification and classification 
of the various types of solutions. On this issue, some progress has been made, particularly for gap 
solitons that bifurcate from linear Bloch modes at the edges of Bloch bands. Specifically, in ID, 
only two gap-soliton families bifurcate from each edge of a Bloch band under self-focusing or self- 
defocusing nonlinearity (Neshev et al. 2003; Pelinovsky et al. 2004; Hwang et al. 2011), and 
the positions of these solitons relative to the periodic medium (or potential) are determined by a 
certain recurrence relation (Hwang et al. 2011). In 2D, four gap-soliton families (or a multiple 
of four families) bifurcate from each edge of a Bloch band under self-focusing or self-defocusing 
nonlinearity (Yang 2010). 

However, for the majority of gap solitons, that do not bifurcate from band edges, systematic 
theoretical treatment is lacking at present. Numerical results indicate that, typically, those so- 
lution families feature multiple branches which bifurcate near band edges at small (but finite) 
amplitude (see Yang 2010 and references therein), yet there has been no analytical explanation for 
this phenomenon. 

In this article, we employ the NLS equation with a spatially periodic potential to analytically 
construct and classify a broad class of gap solitons that bifurcate away from band edges in ID 
periodic media. Our approach is based on the observation that the two families of fundamental 
gap solitons that bifurcate from a band edge can be located only at specific positions relative 
to the underlying potential; otherwise, they would be nonlocal owing to the presence of growing 
tails of exponentially-small amplitude with respect to the peak wave amplitude. However, by 
piecing together such nonlocal solitary waves, it is possible to construct higher-order gap solitons, 
or bound states. The proposed asymptotic theory is illustrated by working out details for bound 
states involving two nonlocal solitary waves. We show that a countable set of bound-state families 
can be constructed. Each family is characterized by the separation distance of the two solitary 
waves involved, and features three distinct solution branches that bifurcate near band edges. The 
analytical predictions are verified by comparison with numerical results. 



2 Preliminaries 

Our study of nonlinear wave phenomena in periodic media is based on the ID NLS equation, 

M t + ®xx - V(x)V + <7# 2 tf* = (1) 

with a periodic potential V(x) and self- focusing (a = 1) or self-defocusing (a = —1) cubic nonlin- 
earity. This equation is the appropriate mathematical model for Bose-Einstein condensates loaded 
in optical lattices (Dalfovo et al. 1999; Morsch & Oberthaler 2006) and laser beam transmission in 
photonic lattices under the paraxial approximation (Yang 2010). Although it is possible to consider 
a general periodic potential V(x) as in Hwang et al. (2011), here, for simplicity, we shall work with 
the sinusoidal potential 

V(x) = V sm 2 x , (2) 

which is 7r-periodic and also symmetric in x, Vq being the potential depth. This potential frequently 
arises in nonlinear optics and Bose-Einstein condensates. 
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Solitary-wave solutions of (1) are sought in the form 



(3) 



where [i is the propagation constant, and the amplitude function ip{x) is real-valued and localized 
in space. Inserting (3) into (1), we find that ip(x) satisfies 



For infinitesimal solutions ip(x), the nonlinear term in the above equation drops out. The resulting 
linear version of (4) is a Mathieu-type equation and, by the Bloch-Floquet theory, its wave spectrum 
features a band-gap structure. Specifically, the linear version of (4) admits two linearly independent 
solutions in the form 



with p(x; fi) being 7r-periodic in x. The wave character of these so-called Bloch modes hinges on 
whether the corresponding wavenumber k is real or complex. By requiring k to be real, one then 
obtains an infinite number of Bloch bands for jjl in which the Bloch modes (5) propagate. These 
propagation bands are separated by gaps, where k turns out to be complex, implying evanescent 
behaviour. 

At the edges of Bloch bands, where the modes (5) switch from propagating to evanescent, two 
Bloch modes in the band collide at either = or = ±1, and a single real- valued Bloch mode 
that is either n- or 2-7r-periodic, arises there. This suggests that edges of Bloch bands are possible 
bifurcation points of solitary- wave solutions of the nonlinear equation (4). These solitary waves 
reside inside band gaps and are the so-called gap solitary waves (or gap solitons). 

The bifurcation of ID gap solitons from band edges was discussed by Pelinovsky et al. (2004) 
using a multiple-scale expansion in powers of the wave amplitude, along with a certain constraint 
obeyed by locally confined solutions of (4) . They identified two families of gap solitons that bifurcate 
from each band edge, namely 'on-site' and 'off-site' solitons, depending on whether the soliton is 
centred at a minimum or maximum of the sinusoidal potential (2). 

More recently, Hwang et al. (2011) re-visited the problem of small-amplitude gap solitons near 
band edges, taking a different approach, that is applicable for a general potential V(x). Rather 
than the constraint utilized in Pelinovsky et al. (2004), Hwang et al. (2011) focused on the 
behaviour of the tails of Bloch-wave packets. These tails are controlled by the coupling of the wave 
envelope to the periodic Bloch mode at the band edge, an effect that lies beyond all orders of the 
usual multiple-scale expansion in powers of the wave amplitude. Hwang et al. (2011) carried this 
expansion beyond all orders using techniques of exponential asymptotics (Yang & Akylas 1997), 
for the case of Bloch-wave packets whose envelope features a single hump. It turns out that the 
tails of such wave packets decay as x — > ±oo, and hence gap solitons arise, only for two specific 
locations of the wave envelope relative to the underlying periodic potential. Both these soliton 
families bifurcate from the linear (infinitesimal) periodic Bloch mode at a band edge, and they 
coincide with the on-site and off-site gap solitons found by Pelinovsky et al. (2004) in the case of 
a symmetric periodic potential. 

In the present work, making use of the asymptotic expressions derived in Hwang et al. (2011) 
for the tails of a single- hump Block- wave packet, we shall construct small-amplitude gap-soliton 
families, in the form of bound states, that comprise two or more such wave packets. In contrast to 
the so-called fundamental gap solitons found in Pelinovsky et al. (2004) and Hwang et al. (2011), 



^xx ~ V(x)4> + ^ + = . 



(4) 



p{x; n) = e tkx p(x; fi) , 



(5) 
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these higher-order soliton families bifurcate at small but finite amplitude close to a band edge, and 
they feature multiple branches that do not connect to fundamental soliton families or band edges. 

In preparation for the ensuing analysis, we now summarize the main results of the multiple- 
scale perturbation procedure for a single-hump Bloch-wave packet (Pelinovsky et al. 2004; Hwang 
et al. 2011). Close to a band edge [i = /Uo> sa y> g a P solitons are expected to be weakly nonlinear 
Bloch-wave packets. The solution to equation (4) is expanded in powers of an amplitude parameter, 
< e< 1, 

tp = eV>o + e 2 V>i + e 3 V>2 + ■ ■ ■ , (6) 



along with 
where r\ = ±1, 



At = Pto + lie , (7) 



= A(X)p(x), (8) 



p(x) = p(x;fio), and X = ex is the 'slow' variable of the envelope function A(X). By imposing 
the appropriate solvability condition at 0(e 3 ), it turns out that A(X) satisfies the stationary NLS 
equation 

DA XX + r)A + aaA 3 = , (9) 

where 

D = \ dVdfc 2 |„ , a = J^p\x) dx j Jo 2 V(*) dx . (10) 
The well-known soliton solution of (9) is 

A{X) = asech X ~ X ° , (11) 

where 

0=^/27^, p=y/\D\, (12) 

and Xq = €Xq denotes the position of the peak of the envelope. This solution exists provided 
Di] < and Da > 0; the first of these conditions requires fi to lie in the interior of the band gap, 
while the second condition can be met in the presence of self-focusing (a = 1) nonlinearity if D > 0, 
or self-defocusing (a = —1) nonlinearity if D < 0. 

It is important to note that the envelope equation (9) is translation-invariant, and Xq is a 
free parameter in the solution (11). As a result, gap solitons of equation (4) could be obtained 
regardless of the position of the envelope (11) relative to the underlying periodic potential. This 
conclusion would seem rather suspicious, though, given that the original amplitude equation (4) is 
not translation-invariant. 

This issue was recognized by Pelinovsky et al. (2004), who pointed out that locally confined 
solutions of (4) must obey the constraint 

/oo 
V'(x)ij 2 (x;xo)dx = 0, (13) 
-00 

which can be readily obtained by multiplying (4) with ip x and integrating with respect to x. Upon 
inserting the perturbation solution (6) and the potential (2) into (13), this condition then restricts 
the peak of the envelope to be at either a minimum (xq = 0) or a maximum (xq = vr/2) of the 
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potential (2), corresponding to on-site or off-site gap solitons, respectively. It is also worth noting 
that the so-called Melnikov function M(xq) in (13), which depends on the shift xo of the envelope 
relative to the potential V(x), is exponentially small with respect to e; hence, the constraint (13) 
for possible locations of gap solitons hinges upon information beyond all orders of the two-scale 
expansion (6). 

3 Nonlocal solitary waves 

Another way of reconciling the two-scale expansion (6) with the fact that gap solitons can be placed 
only at specific locations relative to the potential V(x), is by examining the behaviour of the tails 
of Bloch-wave packets near a band edge. Assuming that they have infinitesimal amplitude, these 
tails are governed by the linear version of equation (4). Hence, for fx inside a band gap and close to 
the edge (iq, using the same notation as in (7) and (10), the asymptotic representation of ip(x) away 
from the solitary-wave core is, generically, a linear combination of two evanescent Bloch modes 



to leading order in \fi — hq\ = e 2 <C 1. As expected, the slow exponential growth/decay of these 
modes is consistent with the behaviour of the tails of the envelope function A(X) according to the 
steady NLS equation (9). 

Now, for the purpose of constructing a solitary- wave solution of (4), the left-hand tail must 
involve only ip+(x), ip{x) ~ eC+ip+(x), so that ip — > as x — > — oo. For generic values of C+, when 
one integrates equation (4) from this left-hand tail to the right-hand side, the right-hand tail would 
involve both ip-{x) and ij) + (x), 



so ip(x) is not locally confined. However, the amplitude E+ of the growing tail in (15) is exponen- 
tially small with respect to e (see below), and hence this tail cannot be captured by an expansion in 
powers of e, such as (6). Only when C + takes certain special values can the growing-tail amplitude 
E + vanish, resulting in a true solitary wave solution. This restriction on C + is consistent with the 
integral constraint (13) utilized by Pelinovsky et al. (2004). 

Hwang et al. (2011) computed the amplitude E + by carrying expansion (6) beyond all orders 
of e for a general periodic potential V(x), utilizing an exponential-asymptotics procedure in the 
wavenumber domain (Akylas & Yang 1995; Yang & Akylas 1997). According to this revised 
perturbation analysis, E + , which indeed is exponentially small with respect to e, does vanish for 
two specific positions of the envelope (11) relative to the underlying potential, thus furnishing two 
families of gap solitons that bifurcate at the band edge. For the symmetric periodic potential 
(2), in particular, E + vanishes when the peak of the envelope is placed at xq = or xq = tt/2, 
corresponding to the on-site and off-site gap solitons, as was obtained earlier by Pelinovsky et al. 



As the details of the exponential-asymptotics procedure are rather involved, we shall only quote 
the main results. When the peak of the envelope (11) is not at a minimum or maximum of the 
potential (xq / 0,7r/2), the resulting 'solitary' waves are nonlocal owing to the presence of both 




(14) 



ip(x) ~ eC-ip-(x) + E + ip + (x) 



{x > 1) 



(15) 



(2004). 
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growing and decaying solutions (14) in at least one of the tails. Specifically, assuming that the 
left-hand tail is locally confined in accordance with (6) and (11), 

V> ~ 2eap{x)e e( - x ~ x °^ fS (x -»• -oo) , (16a) 

the right-hand tail features a growing component of exponentially-small amplitude with respect to 
e, in addition to a decaying component analogous to (16a): 

Y> ~ 2eap(x)e- e{x - Xo)/l3 

+ A7rC^e~ nl3/e sin 2x p(x)e eix - x ° )/l3 {x - x > 1/e) . (16b) 

Note that the amplitude of the growing tail in (16b) is proportional to a constant C. As explained 
in Hwang et al. (2011), this constant is determined from solving a certain recurrence relation that 
contains information beyond all orders of expansion (6). In particular, C > for the sinusoidal 
potential (2). As expected, the growing component of the right-hand tail (16b) is absent when 
x = 0,vr/2. 

By utilizing the symmetries of equation (4) as noted below, it is straightforward to deduce from 

(16) the asymptotic behaviour of nonlocal solitary waves whose right-hand tail is locally confined, 

ij) ~ 2eap{x)e-< x - x °W (x -> +oo) , (17a) 
but the left-hand tail comprises a growing and a decaying component: 
i> ~ 2eap(x)e eix - Xo)/l3 

- 47rC^e- 7r/3/e sin 2x p{x)e-< x - Xo)l ^ {x - x < -1/e) . (17b) 

For the symmetric potential (2), equation (4) is invariant with respect to reflection in x {x — > — x), 
and the Bloch wave p(x) at a band edge is either symmetric or antisymmetric in x, p(—x) = ±p(x). 
Since (4) is also invariant with respect to tp — > —ip, the reflected solution ip may thus be written as 

(17) for both symmetric and antisymmetric p(x). The asymptotic expressions (16) and (17) for the 
tails of nonlocal solitary waves are key to constructing new families of locally confined solutions, in 
the form of bound states, as discussed below. 

4 Bound states 

As indicated above, gap solitons in the form of Bloch-wave packets with the 'sech' envelope (11) 
arise only when the peak of the envelope is at a minimum (xq = 0) or a maximum {x$ = ir/2) of 
the potential (2). Apart from these fundamental soliton states, it is possible, however, to construct 
other gap-soliton families by piecing together two or more of the nonlocal solitary waves found in 
§3. The overall strategy for determining these so-called bound states is similar to that followed in 
Yang &; Akylas (1997) for constructing multi-packet solitary- wave solutions of the fifth-order KdV 
equation. Here, we shall work out the details of finding bound states involving two Bloch-wave 
packets. 
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4.1 Matching of tails 



Consider two nonlocal solitary- wave solutions, ip + (x) and ip~(x), whose 'sech'-profile envelope 
functions are centred at Xq and — Xq , respectively, with > 0. In addition, let ^> ± (x) — > as 
x —> ±oo, so the right-hand tail of ip~{x) and the left-hand tail of ip + (x) are nonlocal. Let us define 
the separation between the solitary-wave cores of the two constituent nonlocal waves as 

S = x+ + x Q . (18) 

Then, assuming S > 1/e, (16b) and (17b) are legitimate asymptotic expressions for these tails in 
the 'overlap' region — Xq <C x <C Xq. Specifically, according to (16b), the right-hand tail of ip~(x) 
is given by 

ij)- ~ ±2ea e - £( - x+x ° )//3 p(x) T e~ nl3 ^ sin 2xq e e ( x+x o )IP p [x) , (19) 

and, according to (17b), the left-hand tail of tp + (x) reads 

i) + ~ 2eae e(:r " a; o )/^( x ) - 47r c^e"^ /e sin2x+ e" e(a; "^ )//3 p(x) . (20) 



e 



Note that the upper sign in (19) corresponds to the case where the envelopes of the two solitary 
waves have the same polarity (sign), while the lower sign pertains to the case of opposite polarity. 

Now, in order for these two nonlocal solitary waves to form a bound state, their growing and 
decaying tail components must match smoothly in the overlap region between the two main cores. 
Based on (19) and (20), this requires 

sin2x - = sin2x+ = T ^^ e ^e~^ , (21) 
and, for this matching condition to be met, it is necessary that 



A 



2ttC /?• 



Clearly, owing to the growing exponential e 71 "^ on the left-hand side of (22), this constraint cannot 
be satisfied for e — > when S is finite. As a result, solution families of bound states are expected 
to bifurcate at a finite amplitude e = e c , say, depending on the separation distance S. Below, we 
shall verify this claim and compute e c for the various solution branches. 



4.2 Solution branches 

Attention is now focused on the matching condition (21), in order to delineate the solution branches 
of bound states that bifurcate close to a band edge. 

Consider first the case where the two nonlocal solitary waves have envelopes with the same 
polarity. Then, the upper sign applies in (21), and sin2xg = sin2x^ < 0. This condition can be 
satisfied in two ways: 

(i) x^ = (m + tj)it + d , x^ = (n — ^)tt + d ; (23a) 

(ii) Xq = (m + t;)tt + d , x$ = nir — d, (23b) 
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where m, n are integers and < d < ir/2. Accordingly, the solitary-wave separation (18) corre- 
sponding to these two possibilities is 



(i) S = Nir + 2d; (ii) S = (N + ±)vr , 
where N = m + n. Moreover, the matching condition (21) reads 



sin 2d = 



2irC /? 3 



(24) 



(25) 



with S given by (i) and (ii) in (24). 

It is convenient to label families of bound states by the integer N which, in view of (24), controls 
the separation between the two nonlocal solitary waves that form the bound states. As noted earlier, 
the present theory assumes that these solitary waves are well-separated, namely, S 3> 1/e, which 
in turn requires N ^ 1/e; the validity of this condition will be verified later in this section. 

In preparation for tracing the solution branches of the family of bound states corresponding to 
a given TV, we write 

d = 7t/4 + 5, (26) 
with — 7r/4 < 5 < 7r/4, so that from (24), we have 



(i) S = S + 25; (ii) S = S , 



(27) 



where So = (N + \)tt. The matching condition (25), with S given by (i) and (ii) above, then leads 
to the following two conditions, respectively, 

cos 25 = W(e, N) e~ 2e5/fi , 



where 



cos 25 = W(e,N), 
W(e,N) 



€ e 7T/9/6 e -6So//3 



(28a) 
(28b) 

(29) 



2nC (3 3 

Based on equations (28), given e, one may find the values of 5, and from (27) the corresponding 
separations S, for which bound states are possible. 

Specifically, from equation (28a), it is easy to see that two values of 5 arise for each e above 
a certain threshold, e^, which is the bifurcation point of the bound states obeying (28a). This 
threshold is associated with a double root, 5 = & , of equation (28a) , where 



^[cos25-W(eP,N)e^ )5 ^) 



s=s; 



(i) 



0. 



From this equation and (28a), it follows that 



^^tan- 1 ^//?. 



(30) 



(31) 



Inserting this result into (28a) , then is found by solving 



cos 



,(1) f ,(1) ,(!)<, 

( tan- 1 = , N) exp | - ^ tan" 1 ^ j . 



(32) 
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Figure 1: Schematic diagram for the three branches of the solution curve (5(e). Two branches 
(solid blue) are obtained from solving (28a) and one branch (dashed red) from solving (28b). The 
bifurcation point e£ with 6c > (solid dot) is found by solving (32), and the bifurcation point 

(2) (2) 

e c with 6 C = (open dot) is found by solving (34). (Online version in colour.) 

For N 3> 1, the two solution branches of 5(e), obtained from (28a) for e > Cc , are plotted 
schematically in figure 1. The bifurcation point e c is the turning point of this double-branch 
curve, and 6 C > 0. For e 3> e c j it is clear from (28a) that cos 25 — > 0, and hence the two branches 
approach 5 = ±7r/4, so d = tt/4 + 5 — > n/2,0. Therefore, in this limit, in view of (23a), the 
two solution branches that bifurcate at approach limiting configurations of bound states, with 
both solitary waves being on-site on one branch (as 5 — > 7r/4) and off-site on the other branch (as 
5 — > —it/ '4). In other words, as the solution curve in figure 1 is traced from one branch to the 
other, the two solitary waves making up the bound state transform simultaneously from on-site 
to off-site. In this transition, 5 decreases from 7r/4 to — vr/4, and according to (i) in (27), the 
separation distance S of the two solitary waves decreases from (N + l)ir to Nir. This transition 
behaviour will be verified numerically in §6; see figure 2. 

Turning next to the second possibility of bound states in (27), equation (28b) furnishes two 
values of 5, 

5 = ±^ cos' 1 W(e,N), (33) 

for each e, as long as W(e, N) < 1. Hence, e must exceed a certain threshold, e > e c , where 

ermi 

~2itC P 3 1 



(2) 

W(er,N) = 1. In view of (29), this bifurcation point is determined from the equation 



: In 



(34) 



Here, the two values of 5 in (33) actually correspond to equivalent bound-state configurations (under 
mirror reflection of the x-axis), as can be easily verified from (23b) with d = j-k + 5. Therefore, 
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(28b) describes only one solution branch. For e S> e c , in particular, <5 — > ±7r/4, and this solution 
branch approaches a limiting bound-state configuration that, according to (23b), involves an on-site 
and an off-site solitary wave. 

We remark that the double root 5 = 5^ = of (28b) for e = is also a solution of (28a) 
for the same value of e. Hence, the solution branch (28b) bifurcates from the point 5 = on the 
lower branch of the solutions (28a), which themselves bifurcate at ei \ as illustrated schematically 
in figure 1. As a consequence, the bound-state family corresponding to a certain N comprises 
three distinct solution branches which are connected with each other. The bifurcation points, 

(2) 

and e c , do not coincide; but when N is large, they are close to each other. Indeed, in the limit 
N » 1, it can be readily shown, using (28), (29) and (34), that ~ ~ fiN' 1 / 2 . Thus, for 
large solitary- wave separation {N S> 1), bound-state families bifurcate near a band edge, but never 
exactly at the band edge itself. Moreover, this scaling of and in terms of N is consistent 
with the condition TV 3> 1/e imposed earlier for the purpose of matching the tails (19) and (20) of 
nonlocal solitary waves forming a bound state. 

We now consider the case when the two nonlocal solitary waves participating in a bound state 
have envelopes with opposite polarity, and the lower sign in (21) applies. Hence, sin 2xq = sin 2xq > 
0. This condition can be satisfied in two ways: 

(i) Xq = mir + d , x^ = rnr + d ; (35a) 

(ii) Xq = rmr + d , x^ = (n + ^)ir — d , (35b) 

for some integers m,n and < d < tt/2. As a result, the expressions (24) for the solitary- wave 
separation S = x^ + Xq are also valid here, so the matching condition (21) gives rise to the same 
equations (28) that describe the solution branches 5 = 5(e) of bound states belonging to the family 
labelled by iV = m + n, where d = 7r/4 + (5as before. Specifically, for a given N, the three solution 
branches found earlier arise in this instance as well: two of these branches bifurcate at e = and, 
in the limit e S> represent bound states with both solitary waves being on-site (5 — > — 7r/4) 
or off-site (5 — > vr/4). Marching along this solution curve from the on-site branch to the off-site 
branch, 5 increases from — 7r/4 to 7r/4, and, in view of (35a), the separation distance S of the two 
nonlocal solitary waves increases from Nir to (N + l)ir. (This transition behaviour is also brought 
out by the numerical results presented in §6; see figure 3.) This transition contrasts with the case 
of same envelope polarity examined earlier, where the separation distance S of the two nonlocal 
solitary waves decreases from (N + l)ir to Ntt when marching along the solution curve from the 
on-site branch to the off-site branch. 

(2) 

Finally, according to (28b), (33) and (34), the third solution branch bifurcates at e = e c and, 

(2) 

for e S> e c , represents bound states with one solitary wave on-site and the other off-site. Moreover, 
it follows from (35) that, in the course of continuation of the bound-state solution from the on-site 
branch to this mixed-site branch, one of the two on-site solitary waves remains on-site and does 
not move, while the other on-site solitary wave moves away and becomes off-site. 

5 Power curves 

In applications, gap-soliton solution branches are often described through 

/oo 
ip 2 dx , (36) 
-oo 
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which is commonly referred to as the soliton power. This quantity is also key to understanding the 
stability properties of gap solitons (Vakhitov &; Kolokolov 1973; Shi et al. 2008; Yang 2010). Here, 
we shall derive an analytical expression for the power of small-amplitude bound states near band 
edges, and trace the solution branches found earlier in terms of P. The theoretical predictions will 
be compared against numerical results in §6. 

By virtue of perfect matching of the tails (19) and (20) in the overlap region — Xq <C x <C Xq, 
one may write the following uniformly valid approximation to bound states involving two nonlocal 
solitary waves (to the leading order): 

if} ~ easech{e(x + Xq )/j3}p(x) ± easech{e(x — Xq) / f3}p(x) ; (37) 

here again the upper (lower) sign corresponds to the case where the envelopes of the solitary waves 
have the same (opposite) polarity. Inserting (37) into (36), the power associated with a bound 
state may then be approximated as 

p ~ i + + r ± i , (38) 

where 

/oo 
p 2 (x)sech 2 {e(xTXo)/(3}dx, (39) 
-oo 

/oo 
p 2 (x) sech{e(x + x^ )//?} sech{e(x - x£)//3}dx. (40) 
-oo 

The integrals above can be evaluated by substituting the Fourier series 

oo 

p 2 (x) = q + ^ q m cos 2mx, (41) 

m=l 

for the even 7r-periodic function p 2 (x), into (39) and (40) and then integrating term-by-term. 
Specifically, in the limit e <S 1, eS S> 1, we find 

J ± ~ 2a 2 ef3q + A7ra 2 f3 2 qi cos 2x± e _7r/3/e + • • • , (42) 
/ ~ 8a 2 e 2 gt)5e- 6S ^ + ■ ■ ■ . (43) 

From (25), however, it is clear that e 2 Se~ eS ^ 3> e~ n P/ e ; hence, the second term in (42) is sub- 
dominant in comparison to (43), and (38) yields the following asymptotic expression for P, 

P ~ 4a 2 eq (/3 ± 2eSe~ tS ^) (eS » 1 , e < 1) , (44) 

with 

1 F* 

qo = - P 2 {x)dx. (45) 
n Jo 

Consider first the case of bound states comprising two solitary waves with envelopes of the same 
polarity, where the upper sign in (44) applies. Corresponding to the first possibility for S in (27), 
there are two bound-state solution branches, 5 = (5(e), governed by (28a), and the associated power 
according to (44) is 



P ~ 4a 2 eg 1 /? + 2e(,So + 25)e- e( - So+25)/l3 j , 



(46) 
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In both (46) and (47) , the second term in the brackets derives from the interaction of the tails of 
the nonlocal solitary waves. Moreover, it should be kept in mind that these expressions are valid 
when e< 1 and eS 3> 1; i.e., close to the band edge and for wide separation of the two solitary 
waves forming a bound state. 

Based on the asymptotic formulae (46) and (47), it is possible to deduce the qualitative be- 
haviour of the power curves. We recall that, for bound states of solitary waves with the same 
envelope polarity, 5 decreases from 7r/4 to — 7r/4 as the solution curve (28a) is traversed from the 
on-site to the off-site branch. Since, for given e, the power (46) is a decreasing function of 5 (for 
eS > (3), the power on the off-site branch is then always higher than the power on the on-site 
branch. Similarly, one can see that the power (47) of the mixed-site solution branch (28b) always 
lies between the powers of the on-site and off-site braches. Hence, the power curve associated with 
a bound-state solution family involving solitary waves of the same polarity, is such that the on-site 
branch is always the lower-power branch, the mixed-site branch lies in the middle and the off-site 
branch is the higher-power branch. Also, in the (e, 5) plane, the middle solution branch bifurcates 
from the lower (off-site) branch (see figure 1); thus, on the power curve, the intermediate-power 
(mixed-site) branch bifurcates from the higher-power (off-site) branch. These general features of 
the power curves will be verified numerically in §6 (see figures 2 and 5). 

Next, in the case of bound states with nonlocal solitary waves having envelopes of opposite 
polarity, where the lower sign in (44) applies, a similar calculation based on (35), along with (27), 
yields 



for the single solution branch bifurcating at e c . Based on these power formulae, it is again 
possible to show that, on the associated power curve, the on-site branch is the lower-power branch, 
the mixed-site branch is the immediate-power branch, and the off-site branch is the higher-power 
branch, just as in the case of same envelope polarity above. However, the intermediate-power branch 
now bifurcates from the lower-power branch, which is the opposite of the conclusion reached earlier 
in the same-envelope-polarity case. These theoretical predictions are confirmed by numerical results 
below (see figure 3). 

6 Numerical results 

We now turn to a numerical investigation of bound states in order to make a comparison of numerical 
results against the predictions of the asymptotic theory discussed earlier. To this end, equation (4) 
is solved numerically by the Newton-conjugate- gradient method (Yang 2010), using the sinusoidal 




(48) 



for the two solution branches bifurcating at e c ' , and 




(49) 
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potential (2) with Vq = 6. For this value of the potential depth, when a = 1 (self-focusing 
nonlinearity) , gap solitons bifurcate from the lower edge fj,o = 2.061318 of the first Bloch band, 
where the diffraction coefficient D is positive, and the parameters (12) of the envelope soliton (11) 
take the values 

a = 1.7146, Q = = 0.6594 . (50) 

Moreover, the constant C in the tail expressions (16b) and (17b) is found by solving a certain 
recurrence relation as discussed in Hwang et al. (2011): C = 1.307. On the other hand, when 
a = — 1 (self-defocusing nonlinearity), gap solitons bifurcate from the upper edge fio = 2.266735 of 
the first Bloch band, where the diffraction coefficient is negative, and the soliton parameter values 
are 

a = 1.681, p = ^/\D\ = 0.7669 , (51) 

with C = 3.0133. 

We begin with a few qualitative comparisons between the analysis and numerics. For this 
purpose, we assume self-focusing nonlinearity (a = 1) and consider bound states involving solitary 
waves with envelopes of the same polarity for N = 10. The numerical results for this family of 
bound states are displayed in figure 2. From the power curve in figure 2a, it is seen that this family 
indeed exhibits three distinct solution branches that bifurcate near the band edge, as predicted 
by the theory. On the lower-power branch, the bound state at point c comprises two on-site gap 
solitons which are separated by (N + 1)tt = Utt (see figure 2c); on the higher-power branch, 
the bound state at point e comprises two off-site gap solitons which are separated by Ntt = 10tt 
(see figure 2e). Thus, as one moves from the lower branch to the upper branch along the power 
curve, the two solitary waves forming a bound state transform from on-site to off-site, and their 
separation decreases by one lattice period (tt). On the middle branch, however, the bound state at 
point d comprises an on-site and an off-site gap soliton which are separated by (N + l/2)ir = 10.57T 
(see figure 2d), confirming that this is the mixed-site branch. At the bifurcation point near the 
band edge (tip of the power curve), the bound-state profile is displayed in figure 2f. This bound 
state comprises two low- amplitude Bloch- wave packets which are well separated, as assumed in 
our asymptotic analysis. These features of the power curve and the corresponding bound-state 
profiles are in perfect qualitative agreement with the asymptotic theory. Furthermore, from the 
amplification of the power curve near the bifurcation point, shown in figure 2b, it is clear that 
the middle branch bifurcates from the higher-power branch, again in agreement with the previous 
analysis in the case of bound states with solitary waves of the same envelope polarity. 

Next, we turn to quantitative comparison between the analysis and numerics. For this purpose, 
we choose to consider bound states comprising solitary waves with opposite envelope polarity (the 
results of comparison for bound states with same envelope polarity are very similar). Specifically, 
once again, we assume self-focusing nonlinearity (a = 1) and take N = 10, but the envelopes of the 
two solitary waves in the bound state now have opposite polarity. Figure 3a displays the numerical 
power diagram near the bifurcation point fi c = 2.04442, and figure 36 illustrates the analytical 

(1)2 

power curves given by (48) and (49) near the analytical bifurcation point fi c = fio — e c = 2.04446, 
where = 0.1368 as obtained from equation (32). (To compute the analytical power curves, we 
first solve equations (28) for 5 as a function of e for < e <C 1, and then use e = y^o — //[.) 
Both the numerical and analytical power curves in figure 3a, b exhibit three branches, as expected. 
In addition, the middle (mixed-site) branch now bifurcates from the lower-power (on-site) branch, 
again in agreement with the analysis in the end of §5. Quantitatively, the numerical and analytical 
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Figure 2: (a) Numerical power curve for bound states involving two solitary waves of the same 
envelope polarity with N = 10 and a = 1 (self- focusing nonlinearity) ; the shaded region is the first 
Bloch band. (6) Amplification of (a) near the bifurcation point, (c-f) Soliton profiles at points of 
the power curve marked by the same letters in (a). Shaded stripes represent lattice sites (locations 
of low potentials); stripe separation is equal to the potential period, n. 
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Figure 3: (a,b) Power curves for bound states involving two solitary waves with opposite envelope 
polarity and N = 10, for self-focusing nonlinearity (a = 1): (a) numerical; (b) analytical. The 
analytical solid blue branches are obtained from equations (48) and (28a), while the analytical 
dashed red branch is obtained from equations (49) and (28b). (c,d) Soliton profiles at the points 
(corresponding to /U = 2.0438) of the lower and upper branches marked by the same letters in (a). 
Solid blue curves: numerical; dashed red curves: analytical as given by (37). (Online version in 
colour.) 



power curves are also in reasonable agreement, especially considering that the analytical bifurcation 
point here is £c = 0.1368, which is not all that small; moreover, for this value of e and N = 10, 
eS ~ 4, which is not all that large. Figure 3c, d shows quantitative comparison of numerics and 
analysis for the bound-state profiles on the lower- and higher-power branches at /i = 2.0438 (near 
the bifurcation point). The corresponding analytical bound-state profiles (red dashed curves) are 
given by (37). Here, the locations of the two solitary waves involved in the bound state, namely, x^, 
are determined by (35); these values are dependent on 8(e), which is found by solving equations (28) 
with e = y^/io — fi\. The analytical solution profiles show good agreement with the numerical ones. 

As another quantitative comparison, we examine the dependence of the bifurcation point ec on 
the parameter N that controls the separation distance of the two solitary waves in a bound state. 
The analytical and numerical values of this bifurcation point against iV are displayed in figure 4 for 
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Figure 4: Comparison between analytical and numerical values of the bifurcation point e c at 
various values of N for a = 1 (self- focusing nonlinearity) : (a) same envelope polarity; (6) opposite 
envelope polarity. The analytical e£ is obtained by solving equation (32). 



solitons with the same as well as opposite envelope polarities and self-focusing nonlinearity (a = 1). 
Here, the analytical bifurcation point e£ is obtained from solving equation (32), and this value is 
the same for both cases of envelope polarity (with the same N). For both envelope polarities, the 
numerical values for approach the analytical ones when N — > oo, confirming the asymptotic 
accuracy of our analysis. 

As a final comparison, we assume self-defocusing nonlinearity (cr = —1) and consider bound 
states involving solitary waves with the same envelope polarity. For N = 10, the numerical results 
for this family of bound states are displayed in figure 5. The power curve features three branches 
(figure 5a), and the middle branch bifurcates from the upper branch (figure 5b), in agreement with 
the asymptotic analysis. In addition, the bound state at point c on the lower branch comprises 
two on-site solitons which are separated by (N + l)n = Hit (see figure 5c), and the bound state 
at point d on the upper branch comprises two off-site solitons which are separated by Nir = Wn, 
again consistent with the previous analysis. It may seem puzzling at first sight that the two on-site 
solitons in figure 5c are of opposite signs, given that the envelopes of these solitons ought to have 
the same polarity. To explain this feature, note that at the band edge relevant here, = 2.266735, 
which is the upper edge of the first Bloch band, the Bloch mode p(x) is 27r-periodic, and its adjacent 
peaks have opposite signs, similar to the function cosx (see Yang 2010, §6.1.1). Since the two on- 
site gap solitons in figure 5c are separated by (N + l)ir = ll-zr, the Bloch function p(x) at the 
centres of these solitons then has opposite sign. Thus, even though the envelopes of the two gap 
solitons have the same polarity, the overall gap-soliton profiles, which are products of the Bloch 
function and the envelope function, have opposite sign, as in figure 5c. 



7 Concluding remarks 

In this paper, an asymptotic theory was developed for stationary gap-soliton bound states consisting 
of two fundamental gap solitons in ID periodic media. It was shown that there is a countable set 
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Figure 5: (a) Numerical power curve for bound states involving solitary waves of the same envelope 
polarity with N = 10 and a = — 1 (self-defocusing nonlinearity) ; the shaded region is the first Bloch 
band, (b) Amplification of (a) near the bifurcation point. (c,d) Bound-state profiles at points of 
the power curve marked by the same letters in (a). 



of such bound state families, characterized by the separation distance of the two solitary waves 
involved, and each family features three distinct solution branches that bifurcate near band edges 
at small, but finite, amplitude. Of these three solution branches, one branch contains bound 
states whose fundamental solitons are both on-site; the second branch contains bound states whose 
fundamental solitons are both off-site; and the third branch represents bound states in which 
one fundamental soliton is on-site and the other off-site. The power curves associated with these 
solution branches were computed asymptotically for large solitary-wave separation. On the power 
diagram, the on-site solution branch is always the lower-power branch, the mixed-site branch the 
middle branch, and the off-site branch is the higher-power branch; in addition, the middle branch 
bifurcates from the upper (lower) branch when the envelopes of the two fundamental solitons in 
the bound state have the same (opposite) polarity. These analytical results were compared with 
numerical results and good agreement was obtained. 

There are several common features between the results of this article and those obtained by 
Yang & Akylas (1997) for two-wavepacket bound states in the fifth-order KdV equation. In that 
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case, a countable set of bound states also bifurcate at finite values of the wave amplitude, and 
each solution family also features multiple branches, with asymmetric-wave branches bifurcating 
off symmetric-wave branches. Furthermore, the techniques for constructing bound states in these 
two models follow along the same lines. Even though the current lattice model (4) is not translation 
invariant, while the fifth-order KdV equation is, remarkably, both models admit very similar classes 
of solitary-wave solutions. 

In addition to being of fundamental interest, the bound states discussed here could also prove 
useful in applications as they allow increased flexibility in the profiles of localized nonlinear modes 
within band gaps. In a recent demonstration of image transmission in photonic lattices, in fact, 
Yang et al. (2011) utilized 2D higher-order gap solitons of various shapes, and those gap solitons 
are closely related to the ID bound states constructed in the present paper. For the purpose of 
assessing the potential usefulness of these ID bound states, of course, it is necessary to examine the 
stability properties of the various bound-state families, a question that will be taken up in future 
studies. 

In this paper, attention was focused on bound states consisting of only two fundamental gap soli- 
tons; extending the analysis to more complicated bound states involving three or more fundamental 
solitons, as well as to bound states in 2D periodic media, will be left to future studies. 

Our final remark is that, since the bound states constructed in this paper bifurcate near Bloch- 
band edges at finite amplitude, their power curves always terminate before band edges and never 
reach them (see figures 2a and 5a). This feature of the power curve is shared by some other types 
of gap solitons, such as the truncated-Bloch-wave solitons (Alexander et al. 2006; Yang 2010). 
Whether the present analysis can be applied to those gap solitons needs further investigation. 
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